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Abstract 

Quantum chromo dynamics (QCD) at non-zero isospin chemical potential is studied in a canonical 
approach by analyzing systems of fixed isospin number density. To construct these systems, we 
develop a range of new algorithms for performing the factorially large numbers of Wick contractions 
required in multi-hadron systems. We then use these methods to study systems with the quantum 
numbers of up to 72 7r + 's on three ensembles of gauge configurations with spatial extents L ~ 
2.0, 2.5 and 3.0 fm, and light quark masses corresponding to a pion mass of 390 MeV. The 
ground state energies of these systems are extracted and the volume dependence of these energies 
is utilized to determine the two- and three- body interactions amongst 7r + 's. The systems studied 
correspond to isospin densities of up to pi ~ 9 fm~ 3 and probe isospin chemical potentials, ui, 
in the range m w < fii < 4.5 m n , allowing us to investigate aspects of the QCD phase diagram 
at low temperature and for varying isospin chemical potential. By studying the energy density 
of the system, we provide numerical evidence for the conjectured transition of the system to a 
Bose-Einstein condensed phase at ui > m n . 



I. INTRODUCTION 



An important goal of nuclear physics is to study the interactions and properties of systems 
comprised of large number of hadrons. Nuclear physics is an emergent phenomenon of 
the Standard Model and as this goal requires an understanding of the strong interaction 
dynamics of multi-hadron systems, it necessitates lattice QCD calculations. In recent years, 
preliminary studies of three- and four- baryon systems have been undertaken [TJ |2] and more 
investigations are underway. In addition, systems involving up to twelve ^'s [3111] or twelve 
K ±J s [5] and systems comprised of more than one species [6] have been studied, allowing 
the various two- and three-body interaction parameters of these systems to be determined 
from the energy shift of iV-meson system at finite density. 

The study of systems comprised of large numbers of hadrons can provide vital insight into 
the structure of high density mater, which may exist in the interiors of neutron stars [7], and 
it is also interesting from a purely theoretical point of view to study the rich phase structure 
of QCD. For systems of high isopin density, pi, and non-zero isospin chemical potentials, 
pi, as will be studied here, a complex phase structure has been conjectured [8]. When pi 
reaches the mass of one pion, pions can be produced out of the vacuum and a Bose-Einstein 
condensate (BEC) is expected to form. At asymptotically large pi, the system is known to be 
a colour superconducting BCS-like state and at an intermediate isospin chemical potential, 
a BEC/BCS transition is expected to occur. However, the exact locations and natures of 
these BEC and the BEC/BCS transitions are unknown and can only be determined by 
non-perturbative QCD calculations. QCD studies of systems with finite baryon density are 
hampered by the sign problem resulting from the non-positivity of the determinant of the 
Dirac operator. However for systems with non-zero isospin chemical potential, there is no 
sign problem. By introducing an isospin chemical potential into the QCD action, non-zero 
isospin density systems have been studied in Ref . [9HTT] , showing hints of some aspects of the 
expected phase structure. Non-zero isospin density system can also be studied by the direct 
computation of correlation function of increasing numbers of pions [3H5] and the extension 
of these methods is the subject of the current work. 

Calculating correlation functions involving many-meson systems (here we will focus on 
many tt + systems) involves computing all possible contractions between quark field oper- 
ators, the number of which naively grows as N\N\ for mesonic systems. Even considering 
symmetries between up and down quarks and identifying vanishing and redundant contri- 
bution, the number of remaining contractions grows exponentially with N. In order to 
overcome this problem, much progress has been made in studying many meson systems in 
Ref. [12] by constructing a recursion relation for correlation functions of systems having 
different number of mesons, taking advantage of the fact that many contractions in the cor- 
relation function of an iV-meson system have been partly computed for an (N — l)-meson 
system. Comparing with direct contractions, the recursion relation [12] saves tremendous 
amount of time, and systems having up to 24 7r + 's have been studied in [13] using it. Since 
the Pauli principle excludes putting more than N C N S = 12 quarks (where N c and N s are the 
number of color and spin components of the quark fields, respectively) in the same source 
location, additional sources are required for iV > 12-meson 1 systems. In additional to re- 
quiring more quark propagators, this complicates the recursion relation and increases the 



For simplicity we will refer to "TV-meson" systems. More correctly, we deal with systems in which a 
conserved quantum number (such as isospin) is fixed to N. 
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computational cost of contractions such that the cost for iV = 24 7r + 's is ~ 100 times that 
of 12 7r + 's. Studying systems of 36 7r + 's becomes extremely time consuming even with the 
recursion relation. 

In the current work, we construct new methods to compute correlation function of sys- 
tems containing large numbers of mesons of one species and also for muti-species systems 
by utilizing the fact that the ground state energy is independent of how the 7r +, s are dis- 
tributed among different source locations [13J. The new methods that are presented herein 
significantly speed up the contractions, and enable us to study even higher density systems, 
that are impractical with other methods. Using one of our new approaches, systems com- 
prised of up to 72 7r +, s are studied on four ensembles of anisotropic clover lattices [H] with 
dimensions 16 3 x 128, 20 3 x 128, 24 3 x 128 and 20 3 x 256. This allows us to investigate 
multi-pion interactions and study the phase structure of QCD at non-zero /x/. In this work, 
we are able to probe the QCD phase diagram from fij = m n up to fij ~ 4.5 m w . We provide 
strong evidence for the Bose-Einstein condensation of the system and attempt to investigate 
the BEC/BCS transition at larger ///. 

The layout of this paper is as follows. In Section. [II], we describe new methods to compute 
the correlation functions corresponding to many pion systems, and compare the performance 
and scalability of each method. Details of the lattice ensembles and our computation of 



the relevant correlation functions in momentum space are discussed in Section. Ill By 



applying the most efficient contraction method, correlation function of systems comprised of 
up to 72 7r + 's are computed, and the ground state energies are extracted in Section. IVl In 



Section. [Vj the two-body and three-body interaction parameters are studied. In Section. [VT 
the QCD phase diagram at non-zero isospin chemical potential is investigated, and the 
transition to a BEC state is identified. 



II. METHODOLOGY OF MULTI-MESON CONTRACTIONS 

Non-zero isospin density meson systems can be studied by evaluating correlation functions 
of many 7r + 's at finite volume (as we work in the context of a relativistic field theory, the pion 
number is ill-defined, however the net isospin of the system is specified in the correlation 
functions below). A correlation function for a system of n = YliLi n i tt + 's with rii 7r + 's from 
source locations (yj, 0) is defined as: 




where the interpolating operator 7r + (x, t) = <i(x, t)7 5 w(x, t) and 7r _ (x, t) = m(x, t)7 5 ci(x, t). 
The correlator C nii ... ;njv can be identified as the term with prefactor YliLi fr° m the 
expansion of det[l + \\Pi + A2-P2 + • • • + A at -Pat], where N is the number of sources, and the 
12 N x 12 N matrices are given by: 
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with 12 x 12 sub-blocks 



Pk,i{t) = S ^ * 0)5 t (x, £; y fc , 0), 



(3) 



where S(x, £;y,0) is a quark propagator between two points. Each is an uncontracted 
correlator describing a quark propagating from source i to source k through the sink at x 
with the quantum number of a 7r + . 

As shown in Ref. [12], a recursion relation for the C nir _ )TlN (t) can be derived by studying 
the properties of the expansion of the above determinant. The C nii ... injv (t)'s have the same 
energy spectrum for all combinations of n«'s as long as n is fixed, so separately computing 
correlation functions of all possible combinations of n^'s is redundant. We can thus iden- 
tify a combined correlator CW^i) as the term having prefactor A™ from the expansion of 
det[l + XA], with 



A = P 1 + P 2 + ... + P 1 
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(4) 



CnT,{t) is a complicated summation of all possible C nijn2i ... jnjv (t) with fixed n, in which we 
do not identify which pions originate at which source. For multiple source contractions, even 
terms representing more than 12 7r +, s located in a single source are included, however such 
terms vanish identically and so do not produce additional noise in numerical calculations. As 
fewer correlation functions are needed, computing Cn n (t) is a computationally simpler task 
than recursively computing all C ni;Tl2r __ ;nN (t). In the following subsections, we will construct 
four algorithms to further speed up the calculation of C n7r (t) and compare each algorithm 
in terms of precision requirement and numerical cost. 



A. Vandermonde Matrix method (VMm) 



As described above, a correlation function of an n-7r + system (C„ 
the coefficient of X n from the power series expansion of det[l + A A] 



det[l + XA] = 1 + \C l7T + A 2 CW + • • • + X 12N C_ 



12Nt, 



can be identified as 



(5) 



where A is a 12 A x 12 A matrix constructed from uncontracted correlators following Eq. (|4]). 



A simple way to get C n7V is by computing Eq. ^ for 12 A different choices of A (Ai, 
The resulting system of equations can be written in the following matrix form 
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The matrix on the RHS of Eq. ^ is a 12A x 12A Vandermonde matrix, for which there 
exist analytical forms for the determinant and inverse (see for example Ref. |15j). The 



4 



inverse matrix then allows us to determine the C n7r 's from the numerical calculation of the 
determinant vector. However, when the number of sources becomes large, elements of this 
matrix can become very small or very large because of the factors of \j> 2 >---> 12N - 1 ^ making 
the computation of the inverse very demanding in precision and eventually resulting in 
significant numerical errors. 

B. FFT method (FFTm) 

By choosing A = exp(z27r/ ■ t) in Eq. (§, the expansion becomes 

det[l + XA] = 1 + e 2i ^°- T C l7T + e^ T C^ + ... + e 24i7rJv ^ T C 12JV7r , (7) 

which contains contributions from signals of frequencies kf , k = 1, 2, . . . 12N. Because of 
this feature, the magnitude of each frequency component can easily be extracted using a 
Fast Fourier Transform (FFT). The magnitude corresponding to frequency kfo is equivalent 
to Ckw times a normalization constant. In order to get better signals, data from multiple r's 
are beneficial, which results in the need to calculate many determinants, in general making 
this method expensive. On the other hand, specific choices of fo and r can minimize the 
number of required determinants. We set r n = n dt, for n = 1,2, ... ,T where dt is the 
minimal time step and T is the closest prime number larger than 12N, and fo = -^-^ and 
then compute det[l + X n A] with A n = exp(z27r/ • r n ). After applying the FFT to this series, 
the amplitude of the frequency kfo is TC^. With such choices of f , r n and T, the number 
of determinants needed to compute is the same as the Improved Combination method (ICm) 
discussed below. 

C. Combination method (Cm) 

The FFTm discussed above is constructed from a certain choice of A's so that the ex- 
pansion of the determinant can be recognized as contributions from different frequencies. 
Similarly, by studying the properties of Eq. (|5]), another choice of A's can be utilized to 
eventually separate det[l + XA] into groups of functions individually depending only on 3 
correlation functions. This method requires us to determine the inverse of a 3 x 3 matrix, 
rather than a 12iV x 12iV Vandermande matrix, to solve for the individual correlators. This 
method is applied by the following steps: 

Step 1: Choose f\ = 1 and compute 

D ( i L \f 1 X) = det[l + f 1 XA]-l. (8) 

Notice that D^(fxX) depends on all correlators Ci n , C 27r , . . . , Cun-k- 

Step 2: Choose / 2 = exp(i7r), and construct the following contractions of the functions 
Di(f n X) to generate the following two new quantities: 

D?\X) = D[ 1 \f 1 X) + f 1 D[ 1 \f2\), 

D?\X) = J Di 1) (/ 1 A) + / 2J Di 1) (/ 2 A). (9) 

By inserting the values of fx = 1, / 2 = —1, it is clear that the D\ (A) only depend on 
C(i+i)7n C(3+i)7r, • • and so the correlation functions have been separated into two groups. 
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Step 3: Choose f 3 = exp(i|), and construct the following combinations of the functions 
Df\f n \) and (f n X): 

D?(X) = D?(X) + f 1 D?(f 3 X), (10) 
£>< S) (A) = D?(X) + f 2 D?(f 3 X), 
D<?\\) = D?\x) + hf 3 Df\f 3 X), 
D<?\\) = D?\X) + f 2 f 3 D?\f 3 X), 

and we see that the D\ (A) for i — 1,2 depends on C(o +2 i) n , C(4+2i)n, ■ ■ ■, and D\ (A) for 
% = 3, 4 depends on C(9_ 2i ) 7r , C(i3_ 2 j) 7r , .... In each step, one function depending on a block 
of CfcTr's is separated into two functions each depending only on half of the C^'s from the 
previous function. We iterate this procedure until blocks of only 3 C^'s are reached. 

To summarize this method, in "step n", f n = exp^^Lj) is chosen, and after this step 

D^ n 1 ^(A), % = l,...,2 n ~ 2 , will be separated into 2 n_1 functions, D^ n \x), each depending 
on 12N/2 n ~ 1 C^'s. Assume Dm~ l \X) is a function depending on a block of C^'s. Two 
functions, -D^m-i anc ^ ^2mi each depending on a half of the original block of C^'s are con- 
structed from Dm ^(A) + g2m-i • ^ _1) (/„-A) and £>£r 1} (A) + 92m • -Dm l \f n • A), where 
the q k s, k = 1,2, ... , 2 n_1 , are prefactors used to construct new functions depending only 
on half of the C^'s, which Dm ^(A) depends on. The prefactor q k in step n is constructed 
in the following way. 
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2 n - 2 + l,2"- 2 + 2,...,2 n - 1 , (11) 

where "Group m" contains 2 m ~ 2 functions for m — 2, 3, . . . , n. This process is repeated until 
functions, D^\\), each depending only on 3 CW's are reached. Eventually det[l + XA] is 
separated into functions, D k n \\), depending on following blocks (B k ): 

Group 1: B 1 — [C^^, Csn-k, C^n-k] 

Group 2: B 2 = [C 2 N-ki CqNtt, CiOATtt] = Csub(Bi)-2iV 

Group 3" | -^3 = [CWtt, CVtVtt, CllAfyr] = Csub(Bi)-JV 

[ B4 — [Cntv, C5N1T, Cww] = Csub(B 2 )-N 

Group n: B k = C Sub(Bfc _ 2 „_ 2) __^, k = 2"~ 2 + 1, 2^ 2 + 2, . . . , 2"" 1 (12) 

where Sub(-Bfc) are the sub indexes of the C's in B k , for example Sub(-Bi) = {4iV, 8N, 12iV} 
and Csub(B!)-2Af = {C 2 Nn, Cqntt, Cwn-k}- The dependence of B k on the corresponding C's 
can be determined from the above recursion relation. 

In order to get the individual C-'s, D k n \\j) is required for three different A/s. Dif- 
ferent choices of A/s have no effect on C i7r 's (we have confirmed this numerically). From 
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FIG. 1: C 2w (t), Ci27r(i) and C23 n (t) calculated from 2 sources by ICm with 64-decimal digital preci- 
sion, denoted as ICm64, and Cm with 64(100)-decimal digital precision, denoted as Cm64(Cml00), 
on a single configuration are compared. Correlation functions from CmlOO agree with those from 
ICm64, however for the same precision, the ICm gives more accurate result than Cm. For C2n(t), 
Cm64 fails because of the numerical inaccuracy. 



the D^ 1 (Aj)'s, the C^'s are extracted by solving the following equation, taking the block 
[CW, C 8Nw , C 12Nw ] for example, 

( D?\Xx)\ / Xf N Xf N X{ 2N \ /CW\ 

D? ] (X 2 ) = Xf Xt N Xl 2N • [ C 8N7T . (13) 
\D^(X 3 )J \XrXl N Xf N J \C l2N J 

Inverting this matrix does not suffer from the numerical instabilities seen in the VMm, 
however as 12N becomes large, even computing the inverse of these 3x3 matrices requires 
high precision. Fig. [T] shows a comparison of the correlation functions computed from 2 
sources by applying the Combination method and Improved Combination method to be 
discussed below. For the Combination method at 64 digit precision, Ci n (t), C 27T (t) and 
C^it) show signs of numerical break down at earlier time slices, which goes away at higher 
precision (100 digit), indicating that even calculating the inverse of the 3x3 matrix needs 
high precision to get correct results. 

As constructed, this method is only applicable to a 2 n source problem. In order to solve 
problems having arbitrary number of sources, we extended this to an Improved Combination 
method in the next section. 



D. Improved Combination method (ICm) 

As there are 12iV terms in the expansion of det[l + XA], the Combination method does 
not allow us to determine functions depending on less than 3 C^'s. A similar problem 
appears in the application of the FFT. In order to use FFT, 2 n data points are required. If 
the number of points in a series is not equal to 2 n , points with value zero must be appended 
to the original series to produce a series of length 2 n . Similarly, we can append additional 
Cfe^-'s to the expansion of det[l + XA], as: 

det[l + XA] = 1 + XC l7T + X 2 C 27T + ... + X 12N C 12N7T + X 12N+1 C {12N+l)n + ... + X 2m C 2m jAA) 

where C p7T = for all p > 12N. The power m is chosen such that 2 m_1 < 12iV < 2 m . 
With this new arrangement, exactly the same prescription discussed for the Combination 
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FIG. 2: Correlation functions on a single configuration at t = 20 from 2 sources computed with 
the Improved Combination method using the arprec library [16J at various precisions: 'arprec X' 
denotes that the calculation is done with X-decimal digit precision. The C nvr (20) for n = 1, 2, ... 24 
all agree for the different precision calculations just as they should, except for the calculation from 
16-digit precision. However C n7r (20) for n = 25, 26, . . . , 32 are all machine zero at each precision. 
The disagreement of 16-digit precision indicates higher precision is needed. A similar comparison 
is shown for the single source correlation functions in the insert. 



Method can be applied, but in the last step the -D^ (A) individually depends only on a 
single correlation function. 

A significant advantage of this method compared with the Cm is that no matrix inver- 
sion is required, so it is consequently less demanding in numerical precision, see Fig. [TJ and 
in addition, problems with arbitrary numbers of sources can be solved with this method. 
Correlation functions appended to the series are solved for simultaneously with the other 
C/cr's, providing a numerical check of the validity of this method. In Fig. [2j correlation func- 
tions calculated from 1-source and 2-sources on a single configuration are shown for different 
precision (we use the "arprec" library [16] to perform arbitrary precision calculations). As 
expected, all C p7r 's for p > 12N are indeed numerically equivalent to zero, decreasing expo- 
nentially as the the numerical precision is increased. Since this method is more numerically 
stable than the Combination method, and can also solve problems of arbitrary number of 
sources, it is used in our further studies. 
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E. Generalization to 2 species from N sources 



The methods discussed above can easily be generalized to two species by studying prop- 
erties of the expansion of det[l + X\A + A2-B], where A and B are uncontracted correlation 
functions of two distinct species, for example ir + and p + . We can write 

det[l + XxA + X 2 B] = 1 + X% + X 1 2 T 1 + ... + X k 2 T k + ..., (15) 

where 

Tk(Xi) = X\CQA,kB + ^ ^ J XiCiA,kB + • • • + ^ ^ j A* f k C(M-k)A,kB, (16) 

where M = 12N is the dimension of the matrices A and B, and the correlation functions, 
C m A,nB, are complicated combinations of correlation functions of a system having m-A's and 
n-B's distributed among different sources in all possible ways. 

The Tj(Ai), for j = 0, 1, . . . , M for one Ai can be separated out by applying the methods 
discussed above with different choices of A2's, and then by applying the method again for 
different choices Ai's for all Tj(Ai)'s, the CmAns's can be separated out. This can be further 
generalized to correlators of arbitrary number of species as necessary. 



F. Performance of different methods 

In order to test the accuracy of different methods at a fixed precision, we com- 
pared correlation functions calculated from the VMm (implemented in MATLAB), the 
FFTm (implemented in MATLAB), and the Cm (implemented in C++ using the "arprec" 
high precision library [16]). We first considered a toy model with matrix elements 
A n>m = sin((m — l)(n — 1) + 2) + zcos(2(n — 1)) for 1 and 2 sources in the top half of Fig. [3j 
For this test, the A's used in the VMm and Cm are randomly chosen between —0.25 and 
0.25, however C nn (t) should be independent of these choices. Results from VMm on 1-source 
agree with those from the FFTm and Cm for any set of A. However for 2 sources, the FFTm 
and the Cm give the same results, but the VMm gives inconsistent results and changes with 
different choices of A's, signaling a breakdown of the VMm and the requirement of higher 
precision. Similar tests have also been performed with the matrix elements A n ^ m extracted 
from real quark propagators and the results are shown in the lower half of Fig. [3j In this 
test, the Recursion Relation method (RRm) has also been used to compute the C nvr 's in 
order to validate the new methods. For the N > 1-source calculation no direct comparison 
with the RRm can be made, since the C R7r computed from the new methods are complicated 
combinations of all C niv .. jniV 's with Y^i=i n i = ™- AVe verified however that the energies 
extracted for these correlators with either method, RRm and Cm, are in agreement. In 
contrast to the toy model, for the real A n ^ m , the VMm gives more accurate results than the 
FFTm. However both tests show that the Cm gives the most accurate results for a fixed 
precision. Tests with real propagators on 2-source shows a break down of Cm on and 
C 2 -ki however this breakdown can easily be corrected by working at higher precision. 

The main purpose of constructing these new methods is to expedite the contractions 
required in computing correlation functions for systems comprised large number of mesons. 
The numerical scaling of the Recursion Relation method, Vandermonde Matrix method, 
Combination method and Improved Combination method (the FFT method costs the 
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FIG. 3: The left panels compare 1-source calculations from the VMm, FFTm and Cm, and the 
right panels compare C nn calculated from 2 sources by the three methods. The real propagator is 
taken from one time slice, t = 20. The Recursion Relation method (RRm) |12j is also compared 
with other methods in the lower left plot as a check on the validity of the Cm. For the 2-source 
calculation in the toy model (top right) with VMm, two different sets of A n s have been used, 
denoted as VMml and VMm2. For VMm applied to the real propagator calculations, only one 
choice of A's is shown. "Cm 16 (32)" denotes that calculation is done using Cm with 16(32) decimal 
digit precision. 



same amount of time as the ICm if /o, r and T are chosen appropriately) are com- 
pared in Table |TJ For each method, we determine how many multiplications are required. 
From Ref . [12] , the computational cost of the recursion relation method is proportional to 
12 4 iV 4 exp(2.8(A r — 1)), where N is the number of sources. The VMm requires a calcula- 
tion of 12 N determinants, one inversion of 12V x 12 N matrix and the multiplication of a 
12iV x 12 N matrix and 12iV x 1 vector. The dominant contribution to the computational 
cost of the other two methods comes from calculating a large number of determinants. For 
the Improved Combination method, a step-n calculation requires the computation of 2 n de- 
terminants, while the Combination method requires 3 • 2 n for a step-n calculation. To solve 
an iV-source problem, the Combination method requires log2(4iV) steps for N = 2 m , where 
m is an integer, and the Improved Combination method requires floor(log 2 (12iV)) + 2 steps. 
Taking account of all the determinant calculations that are needed, and the computational 
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FIG. 4: Comparison of the number of multiplications required for each method (RH axis), and the 
corresponding computation time of C nn (t) for n = 1, 2, . . . 12 N on a single time slice, corresponding 
to one application of the specified contraction method in seconds using a single 2.4 GHz Xeon core 
(LH axis). The computational cost of the ICm is taken from the actual running time, and it is 
used to normalize the time scale so that the projected running time of other methods can be read 
out from the LH axis. 



TABLE I: Scaling of different methods in terms of number of multiplications for an N source 
calculation. 





scaling 


RRm 


12 4 N 4 exp(2.8(A - 1)) 


VMm 


(12iV + 2)(12A0 3 


Cm 


3 . 2 log2(4iV)( 1 2 i V) 3 


ICm 


2 floor(log 2 (12A0)+2( 12 jV) 3 



cost of each determinant (~ (12N) 3 using LU decomposition), the numerical cost of each 
method is tabulated in Table. [TJ and compared in Fig. |4| Although the recursion relation 
method significantly reduces the cost of contractions over the original (12N\) 2 scaling, the 
computational cost of the recursion relation method is much larger than other methods, all 
of which scale similarly. Using the ICm, we now turn to our numerical investigations of 
systems of large number of mesons. 
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TABLE II: Details of the four gauge ensembles with the same lattice space a = 0.1227 ± 0.0008 fm 
used in this paper. iV c fg denotes the number of configurations used in the current calculation. In 
the last two columns, N src is the number of source times used on each configuration and JV mom is 
the number of momentum sources used for each source time. 





L 3 x T (a' 1 ) 


L (fm) 


m^L 


m^T 


^cfg 


iVsrc 




Bl 


16 3 x 128 


2.0 


3.9 


8.8 


180 


8 


33 


B2 


20 3 x 128 


2.5 


4.8 


8.8 


51 


8 


19 


B3 


24 3 x 128 


3.0 


5.8 


8.8 


98 


8 


19 


B4 


20 3 x 256 


2.5 


4.8 


17.6 


147 


16 


7 



III. LATTICE DETAILS 

The calculations in this paper are performed on ensembles of anisotropic gauge field 
configurations with clover-improved fermions [18] that have been generated by the Hadron 
Spectrum Collaboration and the Nuclear Physics with Lattice QCD collaboration. The 
gauge action is a tree-level tadpole-improved Symanzik-improved action, and the fermion 
action [THl |2U] isan/ = 2 + l anisotropic clover action [2T] with two levels of stout smear- 
ing [22J with weight p = 0.14 only in spatial directions (see [2] for more details). In order 
to preserve the ultra-locality of the action in the temporal direction, no smearing is per- 
formed in that direction. Furthermore, the tree-level tadpole-improved Symanzik gauge 
action without a 1 x 2 rectangle in the time direction is used. 

Four ensembles of gauge fields are used in this paper with volumes L 3 x T of 16 3 x 128, 
20 3 x 128, 24 3 x 128 and 20 3 x 256, and with a renormalized anisotropy £ = a s j a t = 3.5, where 
a s (a t ) is the spatial (temporal) lattice spacing. The lattice spacing is the same for each 
ensemble and is a s = 0.1227±0.0008 fm [13], which gives spatial extents L ~ 2.0, 2.5, 3.0 fm 
for L = 16, 20, 24 respectively. The same input light quark mass a t mi = —0.0840 and input 
strange quark mass a t m s = —0.0743 are used in generating each ensemble, giving a pion mass 
of m n ~ 390 MeV and a kaon mass of mx ~ 540 MeV. The quantities m n L and m n T, which 
determine the impact of the finite volume and temporal extent, are m n L ~ 3.86,4.82,5.79 
for L = 16, 20, 24 lattices and m^T ~ 8.82, 17.64 for T = 128, 256, respectively. Details of 
the four ensembles are summarized in Table [TTJ 

In our work, a momentum space representation of the contractions is used and Coulomb 
gauge fixed propagators in time-momentum space, which we refer to as "colorwave propa- 
gators", S u /d(p, T] p', 0), are calculated on each configuration 2 . The colorwave propagator is 
defined as 

S u/d (p, t; p', 0) = e~ iP *Su/ d (x, t; p', 0), (17) 
y 

where 

S u/d (x, t; p', 0) = e Jp ' y ^/d(x, t; y, 0), 
y 



2 As we compute gauge invariant quantities, our results are independent of the gauge fixing procedure. 
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is a solution of the Dirac equation: 

D(y, t; x, t)S u/d {^ t; p\ 0) = e* p '^~ . 



x.f 



The colorwave propagator, S u /a(p, r; p', 0) describes a quark propagating from the source 
(p', 0) to the sink (p, r) in time-momentum space. The colorwave propagator and the 
position space propagator are related by a 3 dimensional Fourier transformation as follows: 

S ttAt (p,t;p',0) = J2 e " PX ^ S u/d (x,t;y,0), (18) 

and the conjugate of a propagator is £*(— p, t; —p', 0) = ^S(p', 0; p, ^75. The 75 hermiticity 
of the colorwave propagator follows from the 75 hermiticity of its counterpart in position 
space. 

A correlation function of one pion with momentum p^ can be constructed by projecting 
both sink and source to the same momentum pf as: 

C lv (jp f ,t) = /^e- 4 ( piX - p2X '^(x / ,t) 75 u(x,t) e^e-^- p ^' u(y,0) l5 d(y',0)\ 

\ x,x' yy' / 

= E (e" ipix e ipy 7 5^(x, f; y, 0) 75 e* P2X ' e-*( p - p ^' ( 75 ^(x', t; y', 0) 7i 

x,x',y,y' 

^ 75(e - P1 x e w Su(X) t; y ; o)) ^ e- l ( p - p ^V P2X ' 7 5(75Sl(x\ t; y', 0) 75 ) \ 

\ x <y x',y' / 

75^„(pi, t; p, 0) • 75(75^(-P2,t; P/ - P, 0)75) ) , (19) 



where Pi — P2 = P/- Each choice of {pi,P2} and p satisfying momentum conservation is 
a separate correlation function with distinct creation and annihilation interpolating fields, 
and we have suppressed the dependence of Ci w on p 1 , p 2 and p. During the calculation, we 
held pi, P2 and Pf fixed and summed over all p's for which we have computed colorwave 
propagators (see Table. [Tl]) in order to get more statistics. In the second step, the definition 
of propagator S u /d(x', t; y , 0) and the 75 hermiticity of the propagator is used. The definition 



of the colorwave propagator, Eq. (18), is applied in the last step. 

The correlation functions of a system having n 7r + 's in a single source, with total mo- 
mentum Pf = n pf can be constructed similarly: 

C nn (t,P f ) = ( (^e^ PlX - p2X ')rf(x',t)7 5 n(x,t) 

\x,x' 

j2 e *py e -(p-p/)^( y , o) 75 d(y\ 0) ) ) , (20) 

vy,y' 



where the dependence of C n7r on p 1; p 2 and p has also been suppressed in Eq. (20). 

Because of the Pauli exclusion principle, systems constructed from a single source in 
momentum space can only reach a maximum of 12 7r + 's. In order to put more pions into a 
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system, additional sources are required. Correlation functions of a iV-source system having 
n = YliLi n^^s with total momentum P/ = 5^t=i P/» are §i ven by 



N 



J] J2 ^e'^^uiyj, 0Hd(y>, 0) ) , (21) 




where rii is the number of pions in the i th source, and momentum conservation 



Si=i(Pi ~~ P2) = YTj=iVfj-> must be satisfied in order for the correlation functions to be 



non- vanishing. The contraction methods discussed in the last section apply equally well in 
momentum space and are used in our work. The elements of the counterpart of uncontracted 
correlation functions defined in Eq. ^ are: 

A k ,i (t) = J2 S (Pi. P) Sj ("Pa. P/« - P) > ( 22 ) 
p 

where k, i label the source and sink, and the dependence on p*, p l 2 , p and is suppressed. 

For the T = 128 (256) ensembles, 8 (16) colorwave propagators are generated on each 
configuration located 16 time slices apart to minimize correlations between propagators. For 
ensembles {51, B2, B3, B4}, {180,51,147,98} configurations and {33,19,19,7} momenta 
are used respectively. In order to reduce contamination from thermal states, a temporal 
extent of T = 256 is desirable for systems of large numbers of pions. On the Bl and 
B3 ensembles, the A ± P (antiperiodic ± periodic propagator) method [231425] is applied 
to effectively double the temporal extent. The validity of this method is investigated by 
comparing results from ensemble B4 (20 3 x 256) and with those from ensemble B2 (20 3 x 128) 
with the A±P method and it is found to be sound at the precision we achieve for the systems 
under consideration as discussed below. 



IV. GROUND STATE ENERGIES 

Previous studies of the energies and isospin chemical potentials [6j [13] on ensemble B2 
showed that thermal states contribute significantly to correlation functions and, even for 
Ci2n{t), the ground state does not dominate in any region of Euclidean time. The expected 
form of correlation functions of an n-7r + system with temporal extent T is [6] 

L™J 

C nn (t) = ( U ) A n m Zle-^-™ +E ^ T ' 2 cosh((£ n _ m - E m )(t - T/2)) + . . . , (23) 

m=0 ^ ' 

where A n m = 1 when m = n/2, otherwise A 7 ^ = 2. E m is the ground state energy of a m-7r + 
system, the are the overlap factors for contribution with m 7r's propagating backward 
around the temporal boundary, and the ellipsis denotes contributions from excited states. 
The ground state contribution comes from the m = term, and thermal states are from the 
m^fl terms in the sum, corresponding to contributions where m 7r + 's propagate backwards 
from the source to the sink around the temporal boundary. For the T = 128 B2 ensemble, 
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FIG. 5: The green data is the effective mass calculated from the original data from ensemble B2, 
and blue line is reconstructed from the ground state energies extracted from the ensemble B4 as 
discussed in the main text. The red line is the fitted value of E nn extracted from the correlators 
of ensemble B4. 



effective mass plots are shown in Fig. [5]for various n, and it is clear that correlation functions 
receive significant contributions from thermal states. Their analysis requires a fit including 
all thermal states, Eq. (23), in order to extract the ground state energy. Since the number 
of free parameters in the fit grows with n, the systematic uncertainty of E nn becomes large 
and we are unable to extract any accurate information at large n. In order to minimize 
contributions from thermal states, a longer temporal extent is required. 

Thermal effects are exponentially suppressed by the larger temporal extent and the en- 
semble with T = 256 has greatly reduced contamination, and a simple single exponential fit 
at intermediate times is sufficient to extract ground state energies, even for E 72n , as shown 
in Fig. |6j Effective mass plots of C 2 o n , CW and C 72 w for this ensemble all show a plateau 
region, and a single exponential fit, only including the term in Eq. (23) with m = 0, is 
enough to extract the ground state energy E n7T . However, for significantly larger numbers 
of pions, a still larger temporal extents would again be necessary. 



A. Energies from 20 3 x 256 ensemble 



Correlation functions, defined in Eq. (21), for systems with the quantum numbers of up 
to 72 7r + 's have been computed on the B4 ensemble. In this paper, only systems having zero 
center of mass momentum are investigated. For a discussion of results for different total 
momenta, see Ref [13] . Because of precision issues, we have computed correlation functions 
from 2, 4, and 6 sources, from which E^^in, £25^-5.48^ and E^^^-k have been extracted 
respectively, where E nw is the ground state energy of a n-n + system at rest. Fig. [7] shows 
CWW, 6*407,- (t) and C 707T (t) from 6-source contractions. The breakdown at earlier time slices 
of C 2 o,t(£) indicates that computations with higher precision are required. Computations 
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FIG. 6: The effective mass of C2o n (t) from the 2-source calculation on the ensemble B4 is shown 
on the left along with the extracted ground state energy represented as a black band. Similarly, 
the effective mass of C^it) (C^it)) from the 4 (6) source calculation on the same ensemble and 
the corresponding extracted ground state energy is shown in the middle (on the right). 



with arbitrary precisions are accessible with the "arprec" library [16], however at the same 
precision, they are ~ 5 times more expensive than with the fixed quad-double precision 
(implemented using the "qd" library [17]). In our main studies, we perform all contractions 
in quad-double precision, and multiply the uncontracted propagators by a prefactor before 
performing the contractions such that the particular C nn (t)'s that we focus on do not suffer 
from the limit of the floating point dynamical range of Quad-double precision (this prefactor 
is removed at the end of the calculation). 

As the correlation functions of systems containing many pions span a large numerical 
range, 10 250 ~ 10~ 250 for C-jo^t) for example, inverting the correlation matrix during a 
correlated fit brings in significant instabilities, thus E nn for n = 1, 2, . . . 72 are extracted 
from uncorrelated fits in this study. The fitting window is chosen between time slices where 
a clear plateau region of the effective mass plot can be seen. Statistical uncertainties are 
constructed from fits to multiple bootstrap resamplings of the ensemble (we use N s = 88 
samples), and systematic uncertainties are estimated by shifting the fitting window forward 
and backward two time slices. 

Since the ground state energy of a system containing many pions becomes large, even 
fitting correlation functions with only one exponential becomes problematic because of pre- 
cision. Taking the 25-7r + system for example, the ground state energy of this system is 
-E-257T = 2.76 in temporal lattice units, and the fit is performed between t/a t = [15,58] ± 2. 
The correlation function varies over 140 orders of magnitude from t — 15 to t — 58. Such 
a large change in magnitude requires care with precision and in order to ameliorate this 
problem, instead of fitting correlation functions directly, we fit the following preconditioned 
correlation functions: 

C' nn (t) = Z' n ex P (5E n t)C rm (t), (24) 

where C nn (t) is the original correlation function, and Z' n , and 5E n are fixed numbers, chosen 
so that C ; nv (t) changes less dramatically inside the fitting window. Since the original corre- 
lation function behaves like a single exponential inside the plateau region where the ground 
state dominates, multiplying another exponential will not change this feature. Furthermore, 
the extracted ground state energy should have no dependence on Z' n and 5E n , which is nu- 
merically confirmed. The preconditioned correlation functions and the corresponding single 
exponential fits for n = 20,40 and 72 are shown in Fig. [8j 



16 




FIG. 7: The correlation functions, C20n(t), C^o n (t) and CfQ^it), calculated from 6-sources with 
quad-double precision and double precision are compared in the left, center, and right plots respec- 
tively. The same calculations done with double precision shows even more severely breakdown, 
indicating that high precision is needed in order to study many pion systems. Although C^ott from 
6-sources with quad-double precision breaks down at earlier time slices, the rescaled C^ott from 
2-source computations, which is shown also in the left plot, is free from precision issues and is used 
in extracting the £>207r- 

B. Antiperiodic ± Periodic propagator method (A ± P method) 

By keeping all factors the same as the ground state Zq extracted from the B4 en- 
semble, we have reconstructed the correlators corresponding to the B2 ensemble by utilizing 
the ground state energies extracted from the B4 ensemble 3 . In Fig. [HJ the reconstructed 
effective masses are compared with those from the correlation functions computed from the 
B2 ensemble, showing agreement within uncertainties. The contamination from the thermal 
states on the T = 128 (B2) ensemble can clearly be seen in the rate at which the plateau 
region (where the ground state energy dominates) shrinks as n increases. For systems with 
a large number of pions, excited states have not died out before thermal states become 
important. 

Since a temporal extent T > 128 is required to get a clean signal for many-pion ground 
state energies, we have investigated the use of the A ± P method (combining propagators 
that satisfy anti-periodic and periodic boundary conditions in the temporal direction to 



While we do not expect = Zq for all m because of the effects of pion interactions, deviations are 



expected to be small (This is also supported by thermal fits using Eq. (23) for small n.) 
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FIG. 8: C' 2 Q n (t) is shown on the left, where the green points are data, the red line is constructed from 
the fit, and two vertical dashed lines indicate the fitting window. Similar plots of preconditioned 
C' i0n (t) and Cy 27r (i) are also shown. 
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FIG. 9: Effective mass plots for 24n + and 487r + correlators. The green data are from ensemble B4 
and the red data are from the A ± P method on ensemble B2. Effective mass plots are consistent 
between these two calculations for all n tt + systems. 



cancel certain modes [23H25] ). On the T = 128 B2 ensemble, we check the validity of 
this method in comparison to the B4 (T = 256) ensemble. In order to see the deviation 
of this method compared with those calculated directly from the T = 256 ensemble with 
anti-periodic boundary conditions in the temporal direction, effective mass plots from the 
two ensembles are compared in Fig. [9j and the ratio of correlation functions from these two 
methods are shown in Fig. 10 The A±P method relays on the exact cancellation of thermal 



contributions, and is seen to work very well 1 n + system, see Fig. 10 For systems with more 



than 1 7r + , the A ± P method starts fail at later time slices, however it still gives consistent 
results at earlier time slices, where ground state energies can be extracted. Energies and 
isospin chemical potentials extracted from the A ± P method are compared with those from 



ensemble B4 in Fig. 11, which shows that the disagreement of extracted ground energies 
below 1%, and at our current precision, the A ± P method provides reliable results for the 
correlators we study. This gives us confidence to use the A ± P method for ensembles Bl 
and B3, where we could otherwise not extract ground state energies for large number of 
pions. 
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FIG. 10: The ratio of the correlation function of n tt +, s calculated by using the A ± P method 
on B2 ensemble, C™(i), compared with that from B4 ensemble, C^. 6 (i), for n = 1,3,5,7,11, is 
shown. 



1.010 




FIG. 11: The ground state energies, E nn , extracted from ensemble B2 (E^) with A ± P method 
are compared with those from ensemble B4 (E^) in the left plot, where the ratio of E^f / E^ 
is plotted. The isospin chemical potentials, /i/, at different densities for the two ensembles are 
compared in the right plot. 
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FIG. 12: Effective mass plots with A ± P method on ensemble Bl are shown here. The effective 
mass of C2o-n- it) from the 2-source calculation is shown on the left along with the ground state energy 
represented as a black band. Similarly effective mass plots of from 4-sources and C^nit) from 
6-sources calculations and the extracted ground state energies are shown in the middle and right 
respectively. 




FIG. 13: Effective mass plots with A ± P method on ensemble B3 are shown. 



C. Energies from 16 3 x 128 and 24 3 x 128 ensembles 



As the A ± P method has been validated on the B2 ensemble, systems having up to 
72 7r + 's has also been studied on ensembles Bl and B3 using this method. Effective mass 
plots with extracted ground state energies from ensemble Bl are shown in Fig. 12 and those 
from ensemble B3 are shown in Fig. 13 All calculations are done with the ICm, and ground 
state energies are extracted with the same statistical method as those in the Section IV A 



The extracted ground state energies from all three volumes are shown in Fig. 14 



V. INTERACTION PARAMETERS 

By considering the energy shifts of two particle states in a finite volume, AE = E2—2E1 = 
2a/p 2 + mf — 2 Ei , Liischer derived a relationship between the phase shift, 5(p), and the 
interacting momentum, p = |p|, given by [26, 27] (see also [28]). 



p cot S(p) = 




20 



18 




n(# of pions) 



FIG. 14: The ground state energies of a system of n-7T + (E n7T ) extracted from ensembles Bl (red), 
B3 (green) and B4 (blue) are shown. The black line represents the total energy of n non-interacting 
pions. 



which is valid for momenta below the inelastic threshold. The regulated three-dimensional 
sum, S(x), is 




where the summation is over all triplets of integers j such that | j | < A. 

By performing an expansion in small 1/L, the energy shift of n identical bosons in a finite 
volume, AE n = E n — nE\, has also been studied up to (D(L~ 7 ) in recent work [29H32] . The 
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resulting shift of energies due to both two-body and three-body interactions is given by [30J : 

2 



a 
ttL 



1 + 



a 
ttL 



+ 



a 
a 

7TL 



[I 2 + {2n - 5) J] 

X 3 + (2n - 7)1 J + (5n 2 - 41n + 63) K 
X 4 - 61 2 J + (4 + n - n 2 )./ 2 + 4(27 - 15n + n 2 )X /C 
+ (14n 3 - 227n 2 + 919n - 1043)£ 



+ n C 3 



192 a 5 . 67ra 3 . , _ 

(To + Txn + 77^77 n + 3) X 



Mtt 3 L 7 



M 3 L 7 



(27) 



where m C n = m!/(n!(m — n)!), and the parameter a is the inverse phase shift at the binding 
momentum of the two body system (below we will refer to this as the effective scattering 
length). This is related to the scattering length, a, and the effective range, r, by 



a = a 



2tt 
L3 



a 6 r 1 



7lL 



The geometric constants entering Eq. (27) are: 



X = -8.9136329. 
£ = 6.9458079, 



J = 16.532316, 
% = -4116.2338, 



(28) 



K = 8.4019240, 
Ti = 450.6392. (29) 



-L . 



The three body parameter fj 3 is constructed from the volume dependent but renormal- 
ization group invariant three body interaction parameter, fj^ , the inverse phase shift, a, and 
the effective range, r, as 



V 3 = Vs 1 - 6 



a 

ttL 



727ra 4 r ^ 

J 1 + -ur 1 



where 



r , , 647ra 4 / r- , \ , , rS 96a 4 _ 
t,=vM + -^(zVS-An) log(/iL) - ^<$ms 



(30) 



(31) 



and the renormalized scale dependent coupling is responsible for the three-body in- 

teractions. The renormalization scheme dependent quantity S defined in the Minimal Sub- 
traction scheme is given by 5ms — —185.12506. 



A. Two-body interactions from Liischer's method 

From the energy difference in the 2-7r + system, Ai^ = -EW — 2/77-^, the relative momen- 
tum of each 7r + , p, in the center of mass frame (COM) can be calculated from the dispersion 
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relation. We determine the effective scattering length 4 , 
momenta {Pi}, on each bootstrap ensemble and applying Eq. (25) 



by calculating the interacting 
and we average over 



all ensembles to get the mean value of a, and the statistical uncertainty. The systematic 
uncertainty is determined by averaging the systematic uncertainty of a on each ensemble 
resulting from the systematic uncertainty of the extracted energies from the choice of dif- 
ferent fitting intervals. The extracted effective scattering length for each volume is shown 
Our results are in agreement with the extractions in Ref. [33] from two-body 



in Table III 



systems studied on the same ensembles. 



TABLE III: The effective scattering length (a) from Liischer's method. The first uncertainty is 
statistical uncertainty and the second uncertainty is systematic. 



V 3 x T 


p 2 /ml 


a(fm) 


m^a 


16 3 x 128 


0.0668(45)(1) 


-0.134(7)(5) 


0.263(15)(9) 


20 3 x 256 


0.0301(9) (0) 


-0.122(3)(1) 


0.238(6)(1) 


24 3 x 128 


0.0143(9)(1) 


-0.106(6)(4) 


0.203(12)(7) 


32 3 x 256 a 


0.00678(54)(81) 


-0.114(9)(13) 


0.223(17)(26) 



"Results for this ensemble are taken from Ref. [55] , 



TABLE IV: The effective scattering length (a) and m^f^rj^ extracted from fits to different ranges 
of n. For a fixed ra max , the x 2 /d.o.f. is larger in smaller volumes, indicating that Eq. (27) fails to 
describe systems of high densities. 





n = [3, 5] 


n = [3, 6] 


V 3 x T 


m n a 




X 2 /dof 


ma 




X 2 /dof 


16 3 x 128 


0.260(14)(2) 


0.70(10)(4) 


1.0 


0.261(H)(1) 


0.67(9)(3) 


1.5 


20 3 x 256 


0.234(6)(1) 


0.80(8)(3) 


0.25 


0.235(6)(1) 


0.79(7)(1) 


0.5 


24 3 x 128 


0.209(H)(4) 


1.61(20)(20) 


0.26 


0.209(H)(3) 


1.59(18)(12) 


0.25 




n = [3, 7] 


n = [3, 8] 


V 3 x T 






X 2 /dof 


m n a 


rd—L 


X 2 /dof 


16 3 x 128 


0.262(14)(1) 


0.64(9)(1) 


3.5 


0.263(H)(1) 


0.62(8)(1) 


5.5 


20 3 x 256 


0.235(6)(5) 


0.79(7)(1) 


1.1 


0.235(6)(1) 


0.76(7)(1) 


2.8 


24 3 x 128 


0.211(H)(2) 


1.56(17)(8) 


0.4 


0.210(H)(2) 


1.50(16)(5) 


1.0 



B. Interaction parameters from small a/L expansion 

The dimensionless qualities m n a and m n f^fj^ can be extracted by fitting AE n to the 



large volume expansion of Eq. (27). The fitting strategy is similar to that used in Liischer's 



As discussed above, a is the inverse phase shift at the binding momentum of the two body system, and 



the scattering length in Eq. ( 25 ) uses the Particle Physics sign convention, and it is negative for repulsive 
interactions. 
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FIG. 15: The ma and mf%r] 3 extracted from different fitting windows [n 
fixed and varying n max . 



mm i ^max 



with Tl r 



method by first fitting to each bootstrap ensemble and then computing the distribution of 
fitted parameters in order to get statistical and systematic uncertainties. There are two ways 
to extract ma. One is by fitting only to AE2 using Eq. (27) with the last two lines set to 
zero, and the other way is by fitting mu ltiple AE n 's, with n > 3, and extracting m^f^fj^ at 
the same time as is shown in Table. 



and Fig. 

the later method are chosen from fits with \ 2 ^ 



IV 
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The final a and mf^r] 3 extracted from 



1. We are forced to to use only few body 
systems as the quality of fit rapidly decreases for large numbers of pions. This suggests that 
the weakly interacting pion model of the system that Eq. (27) encodes is becoming less valid. 
Results for the two-body interaction extracted in both ways agree within uncertainties with 
those extracted using Luscher's method, and are shown in Table |V| The original data for 



the AE n ^s and the results from the fits are shown in Fig. 16 



TABLE V: The effective scattering length (a) from small a/L expansion. The symbol "[2]" indicates 



that only AE2 is used in the fitting, and "[3,6]" means that all AE3 to AEq are used. 



V 3 x T 


m n a[2] 


m,ra[3, 6] 


k cot 8/m-n 




16 3 x 128 


0.259(14)(5) 


0.260(14)(2) 


-3.85(21)(3) 


0.70(10)(4) 


20 3 x 256 


0.234(6)(1) 


0.235(6)(5) 


-4.26(11)(10) 


0.79(7)(1) 


24 s x 128 


0.205(12)(5) 


0.210(H)(2) 


-4.78(25)(7) 


1.50(16)(15) 



The effective scattering length, a, extracted from the three volumes depend on the volume. 
With multiple volumes, Eq. (27) can be inverted to extract both the scattering length, a, 
and the effective range, r. In order to do this, we have also used kcot5/m n determined on 
a matching 32 3 x 256 ensemble from Ref. [33J with all lattice parameters the same, except 
for the volume. As the central value of the scattering length from ensemble B3 deviates 
significantly from the trend of the other ensembles (probably due to statistical fluctuations) 
and is lower than the value extracted for this ensemble in Ref. [33J, we exclude this point 
from our fit. We are using the simplest form, k cot 8/m w = —^-^ + ^^(^t); an d neglecting 
higher order shape parameters as our interacting momenta are small. The infinite volume 
results are l/m n a = —4.60(25) and m^r = 22.7(12.3), which agree with the determinations 



of Ref. [33]. The fit is shown in Fig. 17 
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FIG. 16: The energy differences, AE n , are plot as a function of the number of pions, n, where the 
blue points are the original data, the red bands are the fits, and the black bands are the regions 
where the fits are performed. From the left to right, AE n from 16 3 , 20 3 , 24 3 are shown. 
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FIG. 17: The scattering phase shifts from 16 3 , 20 3 , and 24 3 ensembles in this study, are shown 
as the black data points from right to left respectively. The blue data points are the 24 3 and 32 3 
ensemble results from Ref. |33| from right to left respectively. The 24 3 data is excluded in the fit 
as discussed in the main text. The shaded region is the uncertainty and the star is the infinite 
volume result. 



By utilizing the extracted effective range, r, and the effective scattering length, a(L), from 
the three different volumes, from Eq. (30), the volume dependent parameter r\\ , responsible 
for the three-body interactions can be determined for each volume. The extracted values of 
773 are shown in Fig. 18 The dependence of fj 3 on the volume can be rewritten from Eq. (31 ) 
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FIG. 18: The extracted three-body interaction parameter, rj^(L), is plotted as a function of the 
spatial extent of the lattice, L, (black points). The red line shows the expected dependence of fj^ 



on L from Eq. (32) with C = 4.3, which clearly does not provide a good description of the data. 



into a simpler form 



C+ w log(L), 



where C contains contributions independent of L, and a = 647r(3v / 3 — Atx 
We fit rfe to our data to determine C and the best fit is shown in Fig. 



If 



(32) 

-1.48 x 10 3 . 
However the 

X 2 of the fit is poor and it appears that Eq (32) does not effectively explain the volume 
dependence of our data. This might come from the competing of higher order terms O(j^), 
but it also may be a statistical effect. The large value of fj^ for L = 24 is correlated with a 
down shift of the scattering length a. In Ref. |33j, a value of ma = 0.236(18) (27) was found 
for L = 24, which agrees with the value ma = 0.210(16) (5) found above, but with a large 
central value, perhaps indicating a statistical fluctuation. 



VI. QCD PHASE DIAGRAM AT NON-ZERO m 

In Fig. 19 we show the energy density, e = y, determined from the ground state energies, 
E nn that have been computed on each of the three volumes. For a fixed n, the pions are 
forced to be closer to each other in a smaller volume, and the repulsive interactions between 
them become stronger. This drives up the energy of the whole system. The energy densities 
are weakly dependent on the volume, however there are slightly differences as shown in the 
inset of Fig. [19 
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FIG. 19: Energy densities (e) calculated on 3 different volumes are shown as a function of isospin 
density. The blue points are from the 16 3 ensemble, the black ones are from the 20 3 ensemble and 
the pink one are from the 24 3 ensemble. The inset show the slight difference in energy density on 
three ensembles. 

From the extracted ground state energies, the isospin chemical potential can also be 
determined by a backward finite difference, pi (n) = 4^ ~ E n — E n -\. We calculate pi(n) 
on each bootstrap ensemble, which accommodates correlations between E^s extracted on 
the same ensemble, and the systematic uncertainty of the pi(n) from each ensemble is 
evaluated by adding systematic uncertainty from varying the fit ranges used to determine 
E n and E n _i in quadrature. The final systematic uncertainty on pi(n) is from averaging 
the systematic uncertainties of all the bootstrap ensembles, and the statistical uncertainty 
is the standard deviation of the values of pi(n) on the individual bootstrap ensembles. 



In Fig. 20 , the dependence of pi / m n — 1 on the isospin density pi is shown for the three 
volumes. The isospin chemical potential exhibits similar behaviour in all three volumes, 
where they overlap. At small pi, pi increases at an accelerating rate, in agreement with the 
prediction from chiral perturbation theory (%PT) [8J, however at around pj ~ 0.5 fm~ 3 the 
behaviour of pi starts to change, and the accelerating rate gradually decreases, and at even 
higher isospin density the pj starts to flatten off. This change of behaviour of pi indicates 
that the physical state of the system may be altering. 

The expected phase structure of QCD at non-zero isospin chemical potential has been 
discussed in Ref. [8]. At zero temperature, when pj < m n , there is not enough energy to 
excite a pion out of the vacuum. As soon as pi reaches m n , pions can be produced and the 



27 




Piifm 3 ) 

FIG. 20: The isospin chemical potential, fj,j, is plotted as a function of the isospin density, pj, 
from three lattice ensembles, Bl (red), B3 (blue) and B4 (green). The solid black line is from 
expectations of xPT [8. 



r 




FIG. 21: Expected QCD phase diagram following Ref. [8J. Our calculations at a fixed temperature, 
T ~ 20 MeV probe the phase structure along the red dashed line from pj = m n to pi = 4.5 m,. 
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FIG. 22: The ejtsB is plotted as a function of pi/m n . 



system enters a phase with a pion condensate (BEC) via a second order phase transition. 
At asymptotically large values of pi, the attractive nature of one gluon exchange guarantees 
the existence of a colour-superconducting BCS-like state in which quark-anti-quark Cooper 
pairs are formed. At an intermediate value of pi a BEC-BCS crossover is conjectured [8]. 

In this paper, our calculations are performed at a small but nonzero temperature, T ~ 20 
MeV. With the canonical method used in the current calculation, the lowest isospin chemical 
potential that we probe is pi = by definition as we directly add 7r + 's into the system. In 
the smallest volume, for n = 72 7r + 's (the largest value we consider), an isospin density of 
pi ~ 9 fm -3 is achieved, and the phase diagram is explored from pi = up to pi 4.5 
in this paper as shown by the red dashed line in Fig. 21 



In order to investigate the possible phase transition suggested by the behaviour of the 
isospin chemical potential in more detail, we have also compared the extracted energy density 
with the energy density of a cold degenerate system described by a model of weakly inter- 
acting quarks filling their Fermi sphere up to a maximum momentum kp ~ Ep = pr 13^ 
This Stefan-Boltzmann energy density is given by 



(-SB 



N f N r 
~4^ 



-pj 



(33) 



where Nf = 3 and N c = 3. The ratio of e/ess is plotted in Fig. 22, and exhibits similar 
behaviours in all three volumes. The ratio increases from pi = to a peak around pi ~ 
1.3 m n , and then drops and eventually begins to plateau at around pi « 3 m n . Peak posi- 
tions, Pp eak , for each volume identified from Fig. 
for L 

in infinite volume is Pp eak = 1.30(7) m^. The system for pi < 1.3 m n , can be identified as 
a pion gas. When pi ~ Pp ea ki pions start to condense and the system resides in the BEC 



22|are p T peak = {1.20(5), 1.25(5), 1.27(5)} 
{16,20,24} respectively. With an extrapolation linear in 1/L 3 , the peak position 
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state. The plateau beginning to form beyond \ij ~ 3 m^, may indicate a crossover from the 
BEC to BCS state, however higher precision and larger fij is required to make a definite 
statement. Discretization effects also remain to be investigated. 

Two flavour QCD with finite \ij at large T has been investigated in Ref. where a 
finite temperature deconfinement phase transition was identified at fii < m n , however for 
/i/ > TYi-ff no results were presented. In Ref. [TT], the phase diagram of Nf = 4 + 4 QCD was 
investigated at different temperatures and values fij using the grand canonical approach, 
and a phase transition from a pion gas to a BEC state has also been suggested at fij slightly 
higher than m n , in agreement with the results found here. Two color QCD has been studied 
in Ref. |34] , where the authors identified the transition from vacuum to BEC state and the 
BEC/BCS transition. Somewhat interestingly, the ratio of the energy density and its Stefan- 
Boltzmann limit has also been studied (inset of Fig. 1 in Ref. [SI]), showing qualitatively 
similar behaviour to that found in the current study. 

VII. CONCLUSION 

In this work, we have studied lattice QCD at non-zero isospin chemical potential using a 
canonical approach in which we have investigated systems with the quantum numbers of up 
to 72 7r +, s in three lattice volumes, L 3 ~ (2.0, 2.5 and 3.0 fm) 3 at a pion mass of m n ~ 390 
MeV at a single lattice spacing. In order to perform this study, we have developed several 
new methods for performing the requisite Wick contractions of quark field operators. These 
methods are an enormous computational improvement over previous approaches and their 
accuracy and performance have been carefully investigated. 

In our analysis, we have determined the ground state energies of multi-pion systems in 
three different volumes and have used this to extract the isospin chemical potential of the 
states that are produced. In the smallest volumes, systems with isospin chemical potentials 
of up to fij ~ 1600 MeV are created. By considering the energy density as a function of 
the isospin chemical potential, we provide strong evidence for the transition of the system 
from a weakly interacting pion gas to a Bose-Einstein condensed (BEC) phase at fi ~ m n 
as expected from x?T. At higher values of the chemical potential the system is expected to 
transition to a BCS superconducting state and we have sought numerical evidence for this 
but do not have conclusive results. It is interesting to note that the behaviour of the energy 
density as a function of the isospin chemical potential is very similar to that recently found 
in two-colour QCD with a baryon chemical potential by Hands et al. [33] . 

By focusing on few pion systems, we have extracted the two and three pion interactions, 
determining the scattering length, effective range and the renormalisation group invariant 
effective three-body interaction. The scattering parameters were found to be in good agree- 
ment with other recent determinations and we have attempted to investigate the intrinsic 
volume dependence of the renormalisation group invariant three-pion interaction. We have 
also found that as the density increases and the system transitions to a BEC, it can no 
longer be well described in terms of weak few-body interactions. 
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